Free and constrained symplectic integrators for 
numerical general relativity 

Ronny Richter^ and Christian Lubich^ 

■ Mathematisches Institut, Universitat Tiibingen, 

(^ ! Auf der Morgenstelle 10, 72076 Tiibingen, Germany 

O ' ^richter@na.uni-tuebingen.de, ^lubich@na.uni-tuebingen.de 

> ■ November 4, 2008 

We consider symplectic time integrators in numerical General Rela- 
fj ' tivity and discuss both free and constrained evolution schemes. For free 

n^, evolution of ADM-like equations we propose the use of the Stormer-Verlet 

method, a standard symplectic integrator which here is explicit in the com- 
kjq' putationally expensive curvature terms. For the constrained evolution we 

give a formulation of the evolution equations that enforces the momentum 
constraints in a holonomically constrained Hamiltonian system and turns 
the Hamilton constraint function from a weak to a strong invariant of 
the system. This formulation permits the use of the constraint-preserving 
symplectic RATTLE integrator, a constrained version of the Stormer- 
Verlet method. 
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j"*^ . The behavior of the methods is illustrated on two effectively 1+1- 

dimensional versions of Einstein's equations, that allow to investigate a 
perturbed Minkowski problem and the Schwarzschild space-time. We 
compare symplectic and non-symplectic integrators for free evolution, 

p^ ' showing very different numerical behavior for nearly-conserved quanti- 

ties in the perturbed Minkowski problem. Further we compare free and 
constrained evolution, demonstrating in our examples that enforcing the 

k>( I momentum constraints can turn an unstable free evolution into a stable 

, ; ■ constrained evolution. This is demonstrated in the stabilization of a per- 

C^ ' turbed Minkowski problem with Dirac gauge, and in the suppression of 

the propagation of boundary instabilities into the interior of the domain 
in Schwarzschild space-time. 

PACS numbers: 04.25.D-, 04.20.Fy 

1 Introduction 

The Einstein equations of General Relativity have a Hamiltonian formulation 
that arises as a consequence of the Hilbert action principle in a 3-f 1 slicing 
[H [m HOI [331 [5] • The present article deals with numerical methods that respect 
the Hamiltonian structure in the discretization. 



In various areas of scientific computing, such as the dynamics of particle 
accelerators, molecular dynamics, celestial mechanics, quantum dynamics and 
electrodynamics, symplectic integrators for Hamiltonian systems have been es- 
sential for attaining favorable propagation properties (see, e.g., [571 [511 131] and 
references). There have also been a few papers on symplectic integrators for 
free evolution in general relativity [HI M dHl ESI HSj ■ Very recently, in [52] the 
difficulties in devising constraint-preserving symplectic integrators have been 
addressed but not resolved. 

In this paper we present symplectic integrators for free and constrained evo- 
lution in numerical relativity and illustrate and compare their numerical prop- 
erties in numerical experiments with effectively 1-t-l-dimensional versions of 
Einstein's equations. 

In Section 2 we describe the general framework of the Hamiltonian formula- 
tion and its spatial discretization. For the free, ADM-like evolution we propose 
the Stormer-Verlet method in Section 3. This standard symplectic integrator 
(see [Ml [IZl [29]) here is an implicit method, but it is explicit in the compu- 
tationally expensive terms containing the discretized Ricci scalar. Moreover, 
the implicitness is point-wise in space and is readily resolved by simple fixed- 
point iteration. However, in this free evolution scheme the Hamiltonian and 
momentum constraints are not considered and may drift off. 

In Section 4 we give a formulation of the (spatially discretized) Einstein equa- 
tions that enforces the momentum constraints in a holonomically constrained 
Hamiltonian system. This formulation interprets the shift vector as additional 
momentum variables and at the same time fixes it by gauge conditions. In 
the spatially continuous problem, a result in [3] implies that with the enforced 
momentum constraints, the Hamiltonian constraint function satisfies a conser- 
vation law. Even if this property does not extend to the space discretization, it 
is an extra bonus for this momentum-constrained formulation. 

The problem of constraint growth in numerical relativity has also been tack- 
led, using various techniques, in other free and constrained evolution schemes 
(see e.g. [32 1301 [12 HH E])- Stable evolutions of the black hole binary prob- 
lem are possible since 2005 |3I1 ^2\ E [28l HI] • In the approach presented here 
we avoid the growth of the momentum constraints and additionally respect the 
Hamiltonian structure of the equations. 

Although no symplectic integrators are known for general constraints in 
Hamiltonian systems, such integrators do exist for holonomically constrained 
systems. The most basic of these methods is the RATTLE method, a constrained 
version of the Stormer-Verlet method. In Section 5 we discuss its application 
to the holonomically constrained formulation of Section 4. 

Section 6 presents the 1+1-dimensional examples with which we have done 
our numerical experiments: a perturbed Minkowski problem and Schwarzschild 
space-time. We describe the spatial discretization, the Dirac gauge condition, 
and the boundary treatment that we used in our numerical experiments. 

Section 7 compares a non-symplectic Runge-Kutta method and the symplec- 
tic Stormer-Verlet method on a perturbed Minkowski problem, showing very 
different dynamical behavior in the harmonic energies corresponding to different 
modes. These energies are almost-conserved quantities in the time- continuous 
problem and in the symplectic method, but completely fail to be preserved in 
the non-symplectic method. 

Section 8 compares the free and the constrained symplectic schemes on a 



different perturbed Minkowski problem, illustrating that the unstable free evo- 
lution can be stabilized by enforcing the niomcntum constraints. 

Section 9 compares free and constrained symplectic integrators on the Schwarz- 
schild space-time. It turns out that in the constrained scheme, boundary insta- 
bilities do not propagate into the interior of the domain, as opposed to the free 
scheme. 

Our numerical experiments thus illustrate remarkable properties of symplec- 
tic vs. non-symplectic and constrained vs. free integrators that apparently have 
not been addressed in the literature before. 

2 Hamiltonian formulation and space discretiza- 
tion 

2.1 Super-Hamiltonian and constraints 

For classical general relativity a variety of Hamiltonian formulations is known 
(see [2V for a summary of the most popular ones). Our work is based on the 
famous super- Hamiltonian [121 HD] (see also [31 [21] )• It was discovered as a 
preliminary to the quantization of gravity and relies on a 3-1-1 splitting of space- 
time. The geometry is described by the 3-metric of the spatial slices and their 
extrinsic curvature. 

In the super-Hamiltonian, the position variables are provided by the 3-metric 
hij and the corresponding canonical momenta are denoted tt*-' . They are related 
to the extrinsic curvature Kij by tt*-' — \/h(K\h'^^ — K^^y Here h is the deter- 
minant of the metric hij . The super-Hamiltonian takes the form 



n= iaC + P'Ci)d^x (1) 

with freely specifiable functions a and /?*. The densitized lapse a is related to 
the lapse function A^ by a = N/Vh, and the vector with components /3* is the 
shift vector. The functions C and Ci are given by 



C = Tr'H,j--n\nij-hR, (2) 



where R is the Ricci scalar of the metric hij , anqj 

a = -2h,,DkTT'K (3) 

Solutions to the canonical Hamiltonian equations of motion 

are solutions of Einstein's equations if they satisfy the scalar constraint (or 
Hamiltonian constraint) and the vector constraints (or momentum constraints) 

C = , Ci = . (5) 



^Notice that tt-'*' is a tensor density of weight +1. 



If these constraints are satisfied for the initial data, then they remain satisfied 
in the course of the evolution. However, these functions do not satisfy a conser- 
vation law for general initial data. In this sense, C and Ci are weak invariants 
but not strong invariants. In Section U] we will consider a formulation where the 
momentum constraints Ci are enforced (treating /?* as dynamical variables) and 
the Hamilton constraint function C becomes a strong invariant. 

The system @ is the ADM system presented in [3]. It is weakly hyperbolic, 
but not strongly hyperbolic, and hence it has an ill-posed initial value problem. 
For practical applications one will therefore prefer to adopt another Hamiltonian 
formulation of general relativity, like e.g. a generalized harmonic system jl4) . 
However, we focus on the description of applications of symplectic integrators. 
Adopting other formulations will not change our concept, but the details will 
be more complicated. 

For a reasonable comparison of the properties of the different time integra- 
tion methods in the free evolution schemes we will consider an example where 
the unstable modes of the ill-posed system are not excited in the course of the 
simulation (see section [7]). Concerning the constrained evolution scheme (see 
sections [5]) we consider the modified system (|16l) , where the equations (fTi]) as 
well as the vector constraints are explicitly imposed. Unfortunately nothing is 
known yet about the well-posedness of the initial value problem of this con- 
strained system. 

2.2 Discrete Hamiltonian — general form 

To apply numerical methods one will always approximate the continuous func- 
tions hij, tt'^ , /3* and a through objects with finitely many degrees of freedom, 
and introduce discrete derivative operators that act on these finite dimensional 
spaces. With those ingredients a discrete Hamiltonian is derived by replacing 
the continuous functions and derivative operators in H]) with the discrete ones. 
A concrete example of such a finite-difference discretization is given in Sec- 
tion [^21 With some care in defining the discrete canonical momenta, also finite 
element or spectral discretizations yield finite-dimensional canonical Hamilto- 
nian equations of motion. 

We collect the discrete degrees of freedom corresponding to the functions 
hij, TT*-', a, /?* in vectors q, p, a, /3, respectively. The ordering in these vectors 
is chosen such that components corresponding to the same spatial grid point 
are ordered consecutively. 

Since the super-Hamiltonian consists of terms that are either quadratic in 
the momenta tt'^ or linear in both tt*-' and /3* or independent of the momenta, 
any reasonable discretization of the super-Hamiltonian TL will assume the form 
(we ignore the dependence on the discrete densitized lapse ot in the notation) 

^(q, P) = ^P^S(q)p + [/(q) + /3^D(q)p, (6) 

where S(q) and D(q) are matrices of the appropriate dimensions. The matrix 
S(q) is a square and symmetric matrix. In a finite-difference discretization, this 
matrix is block-diagonal with six-dimensional blocks corresponding to the six 
components of tt'^ at a grid point, since the term in the super-Hamiltonian that 
is quadratic in the momenta does not contain spatial derivatives. We then have 



S(q) = blockdiag {S{qe)), where qg contains the six components of the 3-nictric 
at the grid point xi. 

The discrete Hamiltonian is thus quadratic in the momenta, but the func- 
tional dependence on the position variables is more complicated. The compu- 
tationally expensive term with the discretized Ricci scalar is subsumed in the 
potential U{q). 

The canonical equations of motion for this discrete Hamiltonian then read 

q = S(q)p + D(q)^/3 

(7) 
p = -ip^VqS(q)p-VqC/(q)-/3^VqD(q)p. 

Without further ado, the discretized momentum constraints C't = 0, which read 

D(q)p = , (8) 

are not preserved under the evolution of the discretized system ([7]) , nor are the 
discretized Hamilton constraints C — preserved. There may be an exponential 
or even super-exponential drift away from these constraints along solutions of 
the discrete system ([7]). 

3 Free evolution by the Stormer— Verlet method 

A standard symplectic integrator for Hamiltonian systems is the Stormer- Verlet 
scheme, which is particularly prominent in the area of molecular dynamics and 
enjoys a number of remarkable properties; see, e.g., [26J . When applied to 
([7]), a step from old values (q",p") at time i" to values (q"+^,p"+^) at time 
i"+i — t"- + At reads as follows: 

pn+l/2 ^ pn_^(l(p"+l/2)Tv^s(q")p"+V2 + VqC/(q") (9) 

+ /3^VqD(q")p"+l/2' 
+ (D(q"+lf + D(q"f)/3) 

pn+l ^ p"+l/2_^(^i(p«+l/2)T^^g(q„+l)p„+l/2^^^^(q„+l)(^^) 

+ /3^VqD(q"+l)p"+l/2). 

There is only one evaluation per step of the potential U that contains the com- 
putationally expensive gradient of the discretized Ricci scalar. The substeps 
© and (|10p are implicit in p"+i/2 and q""*"^, respectively. They are solved by 
fixed-point iteration, which is local at every grid point. The choice (3 = for 
the shift vector further simplifies the formulas. The last substep pT|) is explicit. 
We will present numerical experiments with this method in Sections [TMl 

The Stormer- Verlet integrator is a second-order method. We remark that 
higher-order symplectic methods are obtained by suitable compositions of steps 
with different step sizes; see, e.g., [371 Section V.3]. 



q"+i = q" + — ((S(q"+l) + S(q"))p"+l/2 (IQ) 



4 A holonomically constrained Hamiltonian for- 
mulation 

Since the constraints are not taken care of in the free evolution ([7]) , there may 
be an uncontrollable drift in the discretized momentum and Hamiltonian con- 
straints. While there exist symplectic integrators for holonomically constrained 
Hamiltonian systems, see [22 Section VH.l] and 29, Chapter 7], there exist no 
general symplectic integrators for systems where the constraints depend on both 
position and momentum variables, as is the case with the momentum and Hamil- 
ton constraints in general relativity. We therefore look for a reformulation of 
the equations of motion for general relativity that enforces the momentum con- 
straints via holonomic constraints and, as an extra benefit, turns the Hamilton 
constraint function from a weak invariant into a strong invariant (i.e., satisfy- 
ing a conservation law). While this reformulation can equally be done on the 
continuous level, we here present it for the discrete equations of motion ([7]) to 
which we will apply a symplectic integrator in the next section. 

We now consider /3 as a dynamical variable and deal with it in two seemingly 
contradictory ways: 

1. We fix /3 by a gauge condition 

g(q) = with invertible matrix A(q) := Vqg(q)'^D(q)'^ . (12) 

Time differentiation of g(q) = and using ^ for q gives 

Vqg(q)^S(q)p + A(q)^ = 0, (13) 

which shows that indeed (3 is determined by the gauge condition, when 
A(q) is invertible. 

A candidate for the choice of the gauge function g is a discretization of 
the Dirac gauge, djihy/^h^^) = 0. With this choice, A(q) is a discretized 
second-order elliptic differential operator. 

2. We consider /3 as new momentum variables canonically conjugate to po- 
sition variables 7 that are not present in the Hamiltonian ([S]), viz., 

H = i/(q, 7; p, /3) = ip^S(q)p + /3^D(q)p + [/(q) , 

which still is quadratic in the momenta (p,/3). 

In the free evolution, the equations of motion for this extended Hamiltonian are 
together with /3 = -V-^H = and 7 = VpH = D(q)p. 

In order to enforce the discrete momentum constraints D(q)p ~ 0, we there- 
fore impose the holonomic constraints (depending only on the position variables 

(q,7)) 

g(q) = 0, 7 = 0. (14) 

The corresponding hidden constraints obtained by time differentiation of these 
constraints and using the unchanged expressions for q ~ S(q)p + D(q)"^/3 (see 
([7])) and 7 = D(q)p then yields 

Vqg(q)^S(q)p + A(q)/3 = 0, D(q)p = . (15) 



This way, the discrete momentum constraints appear as hidden constraints of 
a holonomicaUy constrained Hamihonian system. Moreover, fi is completely 
determined by the first equation in (|15p . 

The full set of equations of motion with Lagrange multipliers A corresponding 
to the holonomic constraints is now 

q = S(q)p + D(q)^/3 

p = -ip^VqS(q)p - VqC/(q) - /3^VqD(q)p - Vqg(q)A 

together with the constraints (fH|) and p3)) and formally also the equations 

7 = D(q)p, /3 = -/x (17) 

with Lagrange multipliers /x corresponding to the constraints 7 = 0. It is to 
this formulation that we will apply a suitable numerical integrator in the next 
section. 

This formulation, which can be similarly given also for the spatially con- 
tinuous problem, does not enforce the Hamilton constraint. However, in the 
continuous case, a result by Anderson & York ^3^ shows that satisfying the 
momentum constraints Ci — implies that the Hamilton constraint function 
satisfies a conservation law 

{dt-Dfi)C = 0, 

where Dfj is the covariant derivative operator in the direction of the shift vector 

13 = (/30. 

5 Constrained evolution by the RATTLE method 

The RATTLE method (0, [571 Section VH.l], [35] Chapter 7]) is an exten- 
sion of the Stormer-Verlet method to holonomicaUy constrained systems. It is 
symplectic and time-reversible, of second order accuracy, and enforces both the 
holonomic and the derived hidden constraints in the numerical solution. When 
applied to pi)l - pT)) . a step of the RATTLE method consists of the following 
equations, which form a nonlinear system for q"+^ and p"+^. 

1. First half-step for the momentum variables: 

„ Af /I 



-1/2 



2 (-(p"+l/2)TVqS(q")p"+l/2 + VqC/(q") (18) 
+ (/3"+l/2)^VqD(q")p"+l/2 + Vqg(q")A"'+" 

^n+l/2 ^ ^"_^^n.+ (19) 



2. Full step for the position variables 

At / 
TV 



q"+l = q" + ^((S(q"+l) + S(q"))p"+l/2 (20) 

n+l\T , Ti(r,'^\T\ Rri+1/2 



+ (D(q"+l)^ + D(q")^)/3 
7"+i = 7"-H^(D(q")+D(q"+i))p"+i/2 (^" = 0) (21) 



3. Position constraints: 

g(q"+i) = (22) 

7"+i = (23) 

4. Second half-step for the momentum variables: 

p"+i = p"+i/2_^(^i(p«+i/2)T^^s(q"+i)p"+i/2 + VqC/(q"+i) (24) 

+ (/3"+l/2)TVqD(q"+l)p"+l/2 + Vqg(q"+l)A" + l-) 
/3"+i = /3«+i/2 _ ^^"+1,- (25) 

5. Momentum constraints: 

Vqg(q"+')^S(q"+l)p"+l+A(q"+l)/3"+l = (26) 

D(q"+l)p"+l = 0. (27) 

Equations ((IH1)-(I131) determine q"+\ and (I111)-(I171) determine p"+^ The equa- 
tions can be solved by an iterative procedure that requires only the solution of 
linear systems with the matrices A(q") and A(q"+^) and their transposes. 

Iterative solution of (|18 p — (|23 p : We start with q" as initial iterate for q""*"^, 
p" for p"+i/2^ ^» for ^"+1/2^ and for A"'+. With these values we first update 
the iterates q"+^ and p"+i/2 through (|18p and ([^U]) respectively, inserting the 
initial iterates at the right-hand sides of these equations. 

Then, for given iterates of q"+\ p"+i/2 and ^"+^/2, equations ^, ^ 
and (|18p yield an equation for A"'^: 

i (D(q") + D(q"+l)) Vqg(q") A"^+ (28) 

= ^(D(q") +D(q"+l)) (p" - ^(i(p"+i/2)Tv^S(q")p"+l/2 



VqC/(q") + (/3"+l/2)^VqD(q")p"+l/2)') 



The matrix on the left-hand side is 0(At)-close to A(q")^ and therefore invert- 
ible under condition (fT^ . 

Equation (P5|) can be interpreted as the requirement to choose A"'"*" such 
that the momentum constraints are satisfied for the next iterate, p^oxt ■ ^^~ 
stead of using (fT5|) to obtain pj^oxt h'om the current iterates we may also 
assume that it is approximately Pnoxt — p""'"^/^ — Ai/2Vqg(q")AA"''^ with a 
small correction AA"'"*" to the Lagrange multipliers. This leads to the following 
equation, 

A(q")^AA"-+ = -^(D(q") +D(q"+l))p"+l/2. (29) 

We solve this approximately for the increment AA"' , and replace the current 
iterate A"^+ := A"'+ + AA"'+. 



Inserting (|20p into ^^ gives an equation for /3"+^". With the current 
iterates in the argument of g, we solve for the increment A/3"^^ " in a step of 
a simpUfied Newton iteration: 

A(q")A/3"+l/2=-i^g(q"+l) (30) 

+ ^Vqg(q"+l)^ (S(q"+l) + S(q")) Vqg(q«)AA"'+ 

and replace (3"+^^^ := ^"+1/2 _^ A/3"+^/2. We then update the iterates of 
pTi+i/2 g^j^^-j q"+i using (fT5|) and ([20]) respectively with the current iterates on 
the right-hand side. We thus have the following schematic iteration cycle: 

A"'+ -^ /3"+i/2 _^ p»+i/2^ q»+i and iterate. 

Solution of (f24l) - (l2711 : Inserting ([H]) into ^7}i yields a linear system for 

A"^^'^ with the matrix A(q"+^)-^. After solving this system we compute p"+^ 
from (|M)l . and then (3"^^ is obtained from solving the linear system (P5|) with 
the matrix A(q"+^). 

We remark that equations (fT^ and (j25p are ignored, since they only de- 
termine approximations to the Lagrange multipliers fi = —f3 that are without 
further interest. 

Symplecticity. The RATTLE method is symplectic with respect to the canon- 
ical symplectic two-form dq A dp + d-y /\ d(3 for the extended phase space of 
(q, 7;p,/3) restricted to the constraint manifold; see ^3 Section VII. 1], [HI 
Chapter 7]. Because of (i7 = 0, the extended symplectic two- form here actually 
reduces to the original canonical symplectic two- form rfq A rfp. 

6 1+1 dimensional test cases 

After the general discussion in the previous sections we now describe in detail 
the problems considered in our numerical experiments. We derive a simplified 
Hamiltonian in Section EUl and in Section f6. 21 we describe the spatial discretiza- 
tion procedure that leads to a discrete Hamiltonian. In Section [Ol we introduce 
the gauge conditions that we use and discuss some problems that are related to 
that topic. We describe a straightforward boundary treatment in Section 16.41 
before we turn to the particular test problems in Section [^31 

6.1 Simplified continuous Hamiltonian 

For a first test we restrict to simple problems. We therefore consider only those 
solutions of Einstein's equations that satisfy the following requirements: 

hij = for i 7^ j, 

h33 = Ch22, (31) 

where either C = 1 or C = sin^ x^ . The latter case corresponds to spherically 
symmetric space-times, and the class of solutions with C = 1 includes a per- 
turbed Minkowski geometry. Both cases will be considered later on. 



It is easy to show that one obtains solutions in cither of these classes if (pij) 
is satisfied at the initial hypersurface Sq, if moreover at Eg 

tt'^ = for i ^ j, 
TT'^{x\x^,x''K-'/^ = n''ix\x^,x'')C^/^ yx\x\x\x^; i,j = 1,2 

^33 ^ ^-1^22^ (32) 

(where C = 1 or C = sin^ x^ in the two cases, respectively) and if the gauge is 
chosen such that everywhere 

f3^{x\x^,x^)=(3^{x\x^,x^) Vx^x^a;^S^ /32 = = /?^ 

a{x^,x'^,x^)C^^^ =a{x\x'^,x^)C^^^ Va;^S^x^x^ (33) 

To summarize, if (155]) is satisfied at every time slice and if (IM]) , (15^ are satisfied 
at So then (PT|) . (l32p are satisfied for all times, as long as Einstein's equations 
hold. 

Because of the last equations in ([3T|) . ([32|l it is natural to define 



h:=l{h22 + C'h33), *:=^22^^^33_ (34) 

It turns out that with these definitions vr is indeed the canonical momentum 
corresponding to h. 

In the spherically symmetric case we now consider the equatorial hypersur- 
face, i.e., x^ = tt/2 and C = 1- Using the new variables it can be shown that 
the following simplified Hamiltonian provides the correct evolution equations 
for the functions /in, vr^^, h and it: 



n= dx^ 



a — TT TT hiihii — TT irhiih 
- a (^dihdih - 2hdfh + hdihdi log(/iii) + 2Chiih 



(35) 



where ^ = 1 in the spherically symmetric case and ^ = if C = 1. 

This Hamiltonian has a similar structure to the super- Hamiltonian ([T]), its 
discretization will therefore also be of the form © . 

6.2 Space discretization in the 1+1 dimensional setting 

Based on (PSJ) we now derive a discrete Hamiltonian. We introduce two uniforrro 
spatial grids {xi, . . . , xjy} and {xi, . . . , a;j\/}|f| staggered such that Xi = {xi + 
Xi+i)/2. We approximate the functions / = a, /in, tt^^, h, tt by piecewise 
constant functions /^, such that 

f^(x) = f, for X, - Aa;/2 < x < x, + Ax/2, 



^The grids are uniform with respect to the spatial coordinate system, i.e. Xi^i — Xi = Ax 
and Xi_|_i — Xi = Ax = Ax. 

^If we apply periodic boundary conditions then M = N, otherwise M = N — 1. 

10 



and then replace the continuous functions in TL by the piecewise constant ones 
[f ^ f^)- For the shift function /3^ the staggered grid {xi, . . . , xm} is the basis 
of this discretization, i.e., 

{I3^)'^{x) = (3, for X, - Ax/2 < x < x, + Aa;/2. 

Discretizing /3^ on the staggered grid leads to better results in the constrained 
evolution. 

Spatial derivatives are approximated by those piecewise constant functions 
whose functional values are obtained through centered finite differencing, i.e.. 



(9i/)^(a:) = (/,+i-/._i)/(2Ax), 



(36) 



for Xi — Aa;/2 < a; < a;^ + Aa;/2. An analogous formula is used also for the 
derivatives of (3^ on the grid {x}. These functions again replace the continuous 
functions dif and dff, respectively, in the continuous Hamiltonian H. 

To get the discretization of the integrand in (|35p one additionally needs to 
take a logarithm and to perform additions as well as multiplications. These 
are pointwise operations where the continuous functions can be replaced by the 
piecewise constant ones easily. Hence, the integrand of the discrete Hamiltonian 
H is again a piecewise constant function. 

Comparing the discrete Hamiltonian for the 1+1 dimensional setting with 
its general form © we identify the potential C/(q) with the discretization of the 
second line in (1351). The matrix S is derived from the kinetic term 



dx^c 



TT TT hiihii — TT nhiih 



(37) 



(38) 



As we discussed in section 12.21 the matrix S is blockdiagonal with each 
block corresponding to a single grid point, because p7p does not contain spatial 
derivatives. Here, in the 1+1-dimensional setting, the blocks are two dimen- 
sional corresponding to the functions tt^^ and tt respectively, and the elements 
of the blocks can be obtained through 



and the matrix D comes from the shift term 

f dx^ (2TT^^hndif]^+Tr^^P^dihn+n(3^dih 



S'^iij^ii = ahiihii, 



"S'ttIItt — —ahiih, 



^TTTT L). 



(39) 



Concerning the matrix D, the integrand in (|38p does contain derivatives. 
Since we use a staggered grid for the discretization of P^ , here D becomes an 

M X 2N matrix of the form 



M ^ N ~1 



M — N (periodic boundary conditions) 






/ 






\ 



(40) 
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where each star represents a row- vector 

D={ Dpi^ii, Dpi^ ). (41) 

6.3 Gauge conditions 

For the constrained evolution scheme discussed in section|4]we additionally need 
to choose a gauge constraint g(q) = 0. In what follows we derive this function 
from the Dirac gauge [20] , a condition that in the continuous problem fixes the 
spatial coordinate system. This gauge condition is formulated in terms of a flat 
background metric, see [TT] for details. 

In the 1+1 dimensional context the Dirac gauge condition becomes 

di (^{x'r^^/'h;^^%^^') = 0, (42) 

where ^ = 1 in the spherically symmetric case and ^ = for the C — 1 class 
of problems. The dependence on x^ for spherically symmetric systems comes 
from the flat background metric, which is the Minkowski metric in spherical 
coordinates there. 

From the continuous gauge condition (jUJ we derive a discrete one as follows. 
First we observe that we obtain equation (j42|l when we require that the variation 
of the following integral with respect to A vanishes 



Tin '■= / dx 



.lAai((xl)-4«/3/^-'/3^2/3y (43) 



This integral is again discretized by the procedure described in Section [ 
that is, the functions /in, h and x^ are approximated by piecewise constant 
functions corresponding to the grid {xi, . . . ,xjv} and A becomes a piecewise 
constant function based on the staggered grid {xi, . . . , xm}- We thus obtain a 
discrete function He of the grid variables q, x and A. Following the procedure 
in the continuous case we take the derivative of He with respect to A as the 
discrete gauge constraints. 

This defines functions g(q), but as discussed in section|4]we additionally need 
that the matrix A(q) is invcrtiblc. It is not a priori clear that this condition 
is satisfied. Therefore we made several spot checks by calculating the singular 
values of A(q). 

It turns out that in the simulations where the computational domain pos- 
sesses boundaries, in particular for spherical symmetry, we indeed find that the 
ratio of the largest to the smallest singular value was at most 2-10^ and the 
ratio of the smallest to the second smallest was always below 5. But with peri- 
odic boundary conditions it turns out that one singular value is more than ten 
orders of magnitude smaller than the others, and hence A(q) must be regarded 
as singular. 

Constraints for periodic boundary conditions. The origin of the non- 
trivial kernel of A(q) can be found already in the continuous 1+1 dimensional 
formulation. There the expression A(q)/3 translates to 



A{h)p^ 







(44) 
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It can be easily checked that with periodic boundary conditions the equations 
dx^A{h)f = 0, I dx^fA{h) {y'^h:^^''^\ = (45) 



hold for any function /. That means the constant functions are in the kernel of 
the transpose operator A(/i)^, and functions of the form ch}/'^hii (c = const.) 
are in the kernel of A(h). If the Dirac gauge is satisfied then h/hu = const., so 
that the constant functions are also in the kernel of A{h). 

This argument can also be extended to other gauge conditions and probably 
to higher dimensions. We do not want to elaborate on this here, we only mention 
that for periodic boundary conditions A{h) will have a non-trivial kernel for any 
pointwise gauge condition of the form g(/iii, h) — 0. 

Now, the singularity of A(q) means that we cannot solve equations ((HI) 
and ([30|) for AA and A^ in the RATTLE scheme. We therefore solve another 
system which enforces the discrete momentum constraints only up to a multiple 
of the spatially constant vector Ijv = (1, • ■ ■ , 1)"^ & ^'^ ■ It turns out that for 
g(q) = this vector is indeed in the kernels of A(q) and A(q)-^, and that the 
system 



A(q) In \ f Af3 \ _ f i 



lAr \ 13 \0 



(46) 



has a unique solution for any right-hand side f . We then solve systems of this 
form instead of (P5|) and (|30p. and take the obtained A A and A/3 as the right 
correction to the Lagrangian multipliers and the shift variables respectively. 

However, with that procedure we do not have control of the mean value 
of the momentum constraints. We see this immediately when we notice that 
l^D(q)(p— Vqg(q)A) = l^D(q)p—l^A(q)'^ A is independent of A (since In 
is in the kernel of A(q)). Hence, the mean value of the momentum constraints 
cannot be brought to zero by this procedure, and instead of satisfying ([5]), we 
can only fulfill PiD(q)p = where Pi = Inxn — IatI^/-^ is the projector to 
the subspace orthogonal to Ijv. 

The full discrete momentum constraints D(q)p = could be enforced by 
adding an extra gauge condition as a discretization of an integral condition 

f c{hii,h)fdx^ =0 

by 

f^c(q) = 0, 
where the function c(q) and the vector f are chosen such that 

l^D(q)VqC(q)f ^ 0. 

It turns out that this condition corresponds to 



Jc{hn,h)difdx^^O 



in the continuous case. We have, however, not included such an extension by 
an integral gauge condition in our numerical experiments. 



13 



ghost zone, ghost zone, 

non dynamical variables computational domain, dynamical variables non dynamical variables 

« X « X II X « X « X « X « X « X « X « X « X « X « X « X « 



• grid {xi, . . . ,a;iv} X grid {xi, . . . , xm} 

Figure 1: The two spatial grids in the computational domain and in the ghost 
zones. 



6.4 Boundary treatment 

We now describe the implementation of boundary conditions. As we focus on 
the investigation of the time evolution, we tried to keep this aspect simple. 

For the tests in the C — 1 class of solutions we impose periodic boundary 
conditions and there is actually no boundary. But in the spherically symmetric 
case it is not possible to impose periodic boundary conditions. 

For our spatial discretization, the time derivatives of the variables at the 
grid point Xi only depend on the variables at the neighboring grid points Xi-2, 
Xi-i, Xi, Xi+i, Xi+2 and Xi-2, Xi-i, Xi, Xi+i. 

Therefore we extend the spatial grids beyond the computational domain, 
that is, we introduce ghost zones to the left and to the right of the boundary 
(see Figure [T]). The function values at the grid points in these ghost zones are 
then treated as non-dynamical variables. In general these function values are 
fixed through the choice of boundary conditions, and here we simply set them 
to values that we read off from a reference solution. 

The discrete Hamiltonian still has the form ©, we only change the interpre- 
tation of the variables. Therefore the equations of motion for the variables in 
the computational domain have the same form as in the case of periodic bound- 
ary conditions and the variables in the ghost zones do not possess evolution 
equations. 

In the general form of the Hamiltonian ^ we then have to distinguish 
between dynamical variables in the computational domain and non-dynamical 
variables in the ghost zones. In particular we must subdivide the matrix D 
such that the non-dynamical character of the ghost variables can be considered. 
Then, D takes the form (cf. equation (j40| ) 



Di 


Di 


& 


D3 


Dint 

D2 


D2 
D4 



D = Di Dint D2 , (47) 



where Dint corresponds to grid points in the computational domain and the 
remaining matrices to grid points in the ghost zones. 

For S the consideration of ghost variables leads to an analogous result, but 
S is block diagonal and hence the only non vanishing "ghost matrices" are S^ 
and S^. These matrices only depend on variables in the ghost zones. Hence, 
they have no influence on the equations of motion and are therefore irrelevant. 
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We collect the ghost matrices in three matrices 



D := ( Di D2 ) 
I D2 



D 



If we denote the variables in the ghost zones as p, q and /3, respectively, 
then the Hamiltonian can be written as 

ff=^P^Si„t(q)p + C/(q,q) 

+ /3^Dint(q, q)p + /3^D(q, q)p + ^^D(q, q)p + ^^D(q, q)p. (49) 

From this Hamiltonian one obtains the equations of motion for q and p, as well 
as the one-step map in the Stormer-Verlet scheme as described in Sections 12.21 
and El respectively. 

We note that with the chosen space discretization, D(q, q) = 0. Moreover, it 
turns out that Dint (q, q) does not depend on q, and so the discrete momentum 
constraints become 

Dint(q)p = 

and thus do not depend on the variables in the ghost zones. This structure is 
due to the discretization of the shift f3^ at the staggered grid {x}. If the shift 
is discretized at {x} then D does not have the bidiagonal structure (001), but 
there are more non-vanishing elements, and the momentum constraints depend 
on q and p as well. 

The grid functions in the ghost zones are specified as 

/(£,i) = /c.t(£,t) (50) 

for each grid function /, each grid point in the ghost zone x and all times i, 
with given reference functions fcxt in the exterior. This does not apply to the 
Lagrange multipliers A, which are defined only on the staggered grid in the 
interior domain. 

In our example of the Schwarzschild space-time given below, the analytically 
known stationary solution is chosen as initial data and exterior data, leaving 
only a numerical dynamics that reflects the stability properties of the various 
numerical schemes. 

6.5 Test scenarios 

In the following sections wc investigate the properties of the Stormer-Verlet and 
the RATTLE method in two different situations. The first one, a perturbed 
Minkowski space-time, is an example for the case C. = 1, and the second one, 
the Schwarzschild space-time, is a spherically symmetric solution. 
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Perturbed Minkowski. The Minkowski metric describes a flat space-time, 
the analytical solution is (with x = x^, y — x^ , z = x'^) 

g = -dt^ + dx^ + dy^ + dz^. (51) 

It is easy to check that for t = const, slicing this solution really is in the class 
([?T|) with C = 1 and thus the numerical schemes we described are applicable. 

Here we perturb the Minkowski initial data. We denote the perturbations e 
and S, and get that at a slice St the independent components of the 3-metric 
are 

hii^l + eu, h = l + d. (52) 

The canonical momenta become 

7r"=5", n^e, (53) 

and without perturbations the slicing density is identical one and the shift vector 
vanishes: 

a = l+e, f3^=e\ (54) 

The spatial grid is chosen as the uniform grid 

""^ ^ ii ' Jn ' ^-^^ * = l,...,Af, (55) 

and we apply periodic boundary conditions. 

Schwarzschild space-time in isotropic coordinates. The Schwarzschild 
space-time in isotropic coordinates is described by the metric (with R — x^, 
9 = x'^, (j) — x^) 

^ = -(MT^i?F + 16i?4 ^'^'^ +Rdn), (56) 

where dQ^ = d9^ + sin^ 6 dcfP'. For t = const, slicing we thus obtain 

{M + 2Rf ~ {M + 2Rf .. 

The extrinsic curvature and the canonical momenta vanish as well as the shift, 
tt" =0, ^ = 0, ^^ = 0, (58) 

and the densitized lapse is 

_ MR'^{2R~M) 



{2R + My 
The spatial grid is chosen as 



(59) 



i?, = 1 + -^^, i^l,...,N. (60) 

To treat the boundaries we continue the uniform spatial grids beyond the bound- 
aries, introducing K grid points in ghost zones to the left and the right. 

The Schwarzschild solution in isotropic coordinates satisfies Dirac gauge. 
This permits us to obtain initial data easily and to compare the numerical 
results with the exact stationary Schwarzschild solution. Applying a numerical 
method with these initial data yields a purely numerical dynamics that reflects 
the stability properties of the various numerical schemes. 
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7 A perturbed Minkowski problem: symplectic 
vs. non-symplectic integrators 

Symplectic Stormer Verlet vs. non-symplectic ICN scheme. Wc com- 
pare the results of the Stormer- Verlet scheme (see Section [3]) with that of a free 
non-symplectic scheme, the iterated Crank Nicholson method (ICN). For a sys- 
tem of ordinary differential equations y = f (i, y), ICN reads 

ki = ^f(i„,y"), 

k2 = ^f(i„ + AV2,y" + ki), 

k3 = Aff(t„ + Ai/2,y" + k2), 
y"+i=y"+k3. (61) 

In both schemes the densitized lapse as well as the shift are chosen not to change 
with time, a — and /3 = 0. 

We apply the schemes to the perturbed Minkowski example of Section 16.51 
where the perturbation is chosen such that the functions denoted s are Gaussian 
functions with width 1/20, center and height 10~^. The functions denoted 6 
are chosen to vanish. That is, we perturb hn, n, a and (3^, but we keep h and 
^11 fj-Qm f^Q unperturbed Minkowski problem. The case where all functions are 
perturbed by random noise, is discussed in the next section. 

Simulation data. The reason to choose ^ = here is the hyperbolicity of the 
considered system. It turns out that the equations of motion that correspond 
to ((55)) are weakly hyperbolic, but not strongly hyperbolic in the sense of [25] . 
However, if we have h = 1 and tt^^ = at some time to then the continuous 
as well as the discrete evolution equations lead to h = I and tt^^ = 0. We can 
hence omit these functions in the analysis of hyperbolicity. The obtained system 
for /ill and tt is then strongly hyperbolic. 

In the simulations we use grids with 51 or 201 grid points. The size of 
the time step is the same as the spatial grid spacing, i.e., At = Ax. The 
computational costs are comparable for both time-stepping methods. 

Simulation results. It turns out that for both integrators the simulations 
are stable at least until t — 1000. We do not see any growing modes, and we 
expect that the simulations remain stable for a much longer time. The discrete 
momentum constraint vanishes identically, because tt^^ as well as dih vanish. 
The discrete Hamilton constraint does not vanish exactly, but it stays almost 
constant at about 5 • 10"^"^ for both schemes. 

To see how the perturbations evolve, we calculate the harmonic energies 
corresponding to hn. Let hn be the vector composed of the grid values of hu, 
and let hii denote the fcth component of its discrete Fourier transform. Let 

further hii be the fcth component of the Fourier transform of hn = "V-^nH. 

Then, the harmonic energy in the fcth mode is Ek — ^(|hiiP -I- fc^|hiip). 

As we see in Figure [2 the E^ are approximately preserved by the Stormer- 
Verlet scheme over long times, whereas with the ICN scheme they incorrectly 
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, N = 51) 



Et (Stormer-Verlet, N = 51) 




Figure 2: Harmonicenergiesin the simulations with 51 grid points. Left: Resuhs 
of the iterated Crank Nicholson scheme, Right: results of the Stormer-Verlet 
scheme. 



decay exponentially until they reach the level of round-off error. The exponential 
decay is much faster for high frequencies. 

Very similar results are obtained in the computations with 201 grid points. 
It turns out that for even numbers of grid points the highest frequency of the 
system grows quadratically for symplectic as well as non-symplectic evolution, 
whereas the remaining Fourier modes show the same behavior as in the case of 
odd numbers of grid points. 

Discussion of results. To sum up the essential feature in this numerical ex- 
periment, we have seen that the harmonic energies are preserved for long times 
by the Stormer-Verlet scheme, whereas they decay quickly in the ICN scheme. 
The symplectic Stormer-Verlet scheme thus reproduces the propagation of per- 
turbations in this problem remarkably better than the ICN scheme or in fact 
any explicit Runge-Kutta scheme. 



8 A perturbed Minkowski problem: stabiliza- 
tion by constraints 

In this section we compare results of the Stormer-Verlet scheme of Section |3] and 
the RATTLE scheme of Section \S\ with constraints imposed as in Section 16.31 
In both schemes we choose the densitized lapse not to change with time, a = 0, 
and in the Stormer-Verlet scheme also the shift is independent of time, /3 = 0. 



Simulation data. Again we apply these schemes to the perturbed Minkowski 
example of Section 16.51 with periodic boundary conditions, but here we choose 
the perturbation functions e as well as S to be noise in the interval (—2.5 • 
10~^/A^^, 2.5-10^^/A^^). In particular, S ^ so that the free evolution equations 
are only weakly hyperbolic. We choose the time step At = Ax, and apply the 
schemes for TV = 50 as well as A^ = 200 grid points. This example is the robust 
stability test suggested in [1] . 
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constraints (N = 50) 
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le+15 r , ,— r- 




Figure 3: Behavior of the numerical solution. Left: The maximum norms of 
the Hamilton and momentum constraints C and Ci, right: the evolution of the 
mean value of hn of hn. 



Simulation results. We first investigate the constraints. We see in Figure 
[3] that the maximum norm of the Hamilton constraint function grows linearly 
in the Stormer-Verlet scheme and is almost constant at about 10^^ for the 
RATTLE scheme. Moreover, in the RATTLE scheme the shape of the Hamilton 
constraint function changes very little with time. 

Concerning the momentum constraint function, in the Stormer-Verlet scheme 
its maximum norm stays constant at about 10^^ until t « 200 and grows su- 
perexponentially afterwards. The RATTLE scheme on the other hand provides 
very small momentum constraints of about lO^^^-lO"^"* that grow cubically 
with time. A closer look reveals that the momentum constraint function in the 
RATTLE scheme is nearly constant already after the first timestep, and the 
cubic growth is due to a growth of the mean value (see Section [6.3p . 

From Figure [3] we also see that in the Stormer-Verlet scheme the error of the 
mean value of hn (we denote this mean value hn) grows with order six in the 
beginning, before a superexponcntial growth becomes important. In the begin- 
ning the errors of the functions hn, h and tt themselves grow cubically, linearly 
and quadratically respectively, while the error of tt^^ stays nearly constant. 
However, at later times the errors of hn and tt^^ grow superexponentially. 

In the RATTLE scheme the error of the mean value hn grows quartically. 
The function hn itself is almost constant in space, and its error is basically due 
to the error of the mean value. It turns out that the errors of the functions h 
and TT grow linearly, and the error of tt^^ stays nearly constant. 

We observe that the main errors of hn in the RATTLE scheme come from 
errors of its mean value. Since it is also interesting to see how the remaining 
components behave, we investigate the deviations from the mean value sepa- 
rately. In Figure [3] we see that the maximal deviation of hn from its mean 
value grows cubically in the beginning for the Stormer-Verlet scheme. Compar- 
ing that deviation in the coarse grid [N = 50) and in the fine grid {N = 200), 
we do not see any difference at all. In the RATTLE scheme the deviation is at 
most 2 • 10"^*^ and it is more than one order of magnitude smaller in the fine 
grid. 

We also look at high-frequency errors of hn. We see in Figure H] that in the 
Stormer-Verlet scheme the component of hn that corresponds to the highest 
frequency grows cubically in time. Moreover it turns out that this component is 
about four times larger in the fine grid than in the coarse grid. In the RATTLE 
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Results of Stormer-Verlet 
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Figure 4: The maximal deviation of hn from its mean value, i.e. ||hii — hii||oo, 
and the highest frequency component of the Fourier transform hn. Left: results 
of the Stormer-Verlet scheme. Right: results of the RATTLE scheme. 



scheme these high frequency errors are very small, and moreover they are about 
four times smaller in the fine grid. 



Discussion of results. We have seen that in the RATTLE schemes the error 
of the Hamilton constraint stays almost constant. This is in agreement with the 
result by Anderson and York [5] that for the continuous problem the Hamilton 
constraint function satisfies a conservation law when the momentum constraints 
are satisfied (cf . also the comment at the end of Section U]) . It is remarkable 
that also in the Stormer-Verlet scheme the Hamilton constraint function grows 
only linearly, although the momentum constraint is not satisfied and grows 
superexponentially at the end of the simulation. 

In the RATTLE scheme we force the momentum constraints to be satisfied, 
but we discussed in section 16.31 that we do not have control about their mean 
value. The numerical results show that even this mean value is small. It is 
interesting that for the RATTLE scheme, the errors of the function hn are also 
mainly due to errors of its mean value. 

In the Stormer-Verlet scheme the high frequency errors grow cubically with 
time and become larger in the fine grid. These errors probably trigger the 
superexponential growth that occur at the end of the simulation. It is typical 
for discretized weakly hyperbolic systems that high frequency errors become 
larger when the grid is refined (see e.g. |16J). High frequency errors do not 
play a significant role in the RATTLE scheme. They remain small and become 
even smaller when more grid points are used for the spatial discretization. This 
example illustrates that imposing constraints may stabilize the problem. 



9 Schwarzschild space-time: free vs. constrained 
scheme 

In this section we again compare results of the Stormer-Verlet and the RATTLE 
method. As in section [5] wc choose the densitized lapse to satisfy <i = in both 
schemes, and in the Stormer-Verlet scheme also /3 = is fulfilled. 
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I (RATTLE) 




Figure 5: The difference between tire numerical and the analytical solution 
h"J"" — hj^^ for the Schwarzschild space-time. Left: results of the Stormer- 
Verlet scheme, right: results of the RATTLE scheme. 



Simulation data. We apply the schemes to the Schwarzschild space-time 
(see Section 1575]) . Again we choose the time step Ai = Ax, and use grids with 
A^ = 51 as well as A^ = 201 grid points. We impose boundary conditions as 
described in Section [ 



Simulation results. From Figure [5] we see that the error of the numerical 
solution grows quadratically in the beginning but at some point before the 
evolution breaks down, a more rapid growth becomes important. The period 
of quadratical growth of errors is for the Stormer-Verlet scheme in the order of 
M (the mass of the black hole) and it is longer for the coarser grid. With the 
RATTLE scheme, the rapid growth of errors starts at about lOM, and moreover 
the errors in the fine grid are always smaller than those in the coarse grid. 

A closer look reveals that in the results of Stormer-Verlet there are high 
frequency errors near the boundary. These errors propagate from the boundary 
into the interior of the computational domain and further amplify. We also 
observed that small perturbations of the initial data as well as the use of, e.g., 
ICN instead of Stormer-Verlet have negligible effects on the solution. In the 
RATTLE scheme high frequency errors of the functions hn, h, tz^^ and tt are 
very small. 

Again we investigate the Hamilton constraint. It turns out that in the 
Stormer-Verlet scheme it behaves analogous to the functions hn, h, tz^^ and tt. 
That is, there are high frequency errors at the boundary that propagate into the 
interior and amplify (see Figure [6|). In Figure [7] we indeed see that the Hamil- 
ton constraint functions in the interior (in the interval R = x^ ^ [1-25, 1.75]) 
and near the boundary are of comparable size. Moreover we see an exponential 
growth of the Hamilton constraint function. 

In the RATTLE scheme we observe a different behavior. There the Hamilton 
constraint function becomes large near the boundaries, too (see Figure [6]), but 
the propagation into the interior is suppressed. If we compare the Hamilton 
constraint function in the interior and in the whole computational domain then 
we see from Figure [7] that it stays almost constant in time in the interior until 
the evolution is on the brink of breaking down. 
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C,N = 201, t = 1.5 (Stormer-Verlet) 

1.5 



|C|, N = 201, f = 1.5 (RATTLE) 
1 




Figure 6: The discrete Hamilton constraint function in the Stornicr-Verlet 
scheme (left) and in the RATTLE scheme (right) at t = 1.5 for iV = 201 
grid points. 
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Figure 7: The maximum of the Hamilton constraint function in an interval 
x^ E I for Schwarzschild space-time and A^ — 201 grid points. We consider 
an interval in the interior of the computational domain, I — [1.25,1.75], and 
the whole computational domain I ~ [1,2]. Left: results of the Stormer-Verlet 
scheme, right: results of the RATTLE scheme. 
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Discussion of results. We have seen that in the free evolution again high 
frequency errors occur, whereas they are absent in the constrained evolution 
scheme. Our simple choice of boundary conditions of Section [Ol leads to bound- 
ary instabilities, and so the main source of errors are the boundaries. In the free 
evolution the errors at the boundaries quickly propagate into the interior of the 
computational domain and amplify. In the constrained scheme, the propagation 
into the interior is suppressed. 
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